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Abstract 

We study theoretically and numerically the steady state diffusion controlled 
reaction A + B ^ ^, where currents J of A and B particles are applied at op- 
posite boundaries. For a reaction rate A, and equal diffusion constants D, we 
find that when AJ"^/^!)^^/^ ^ 1 the reaction front is well described by mean 
field theory. However, for AJ^^/^Z)^^/^ S> 1, the front acquires a Gaussian 
profile - a result of noise induced wandering of the reaction front center. We 
make a theoretical prediction for this profile which is in good agreement with 
simulation. Finally, we investigate the intrinsic (non-wandering) front width 
and find results consistent with scaling and field theoretic predictions. 
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Recently there has been considerable interest in the properties of diffusion limited chemi- 
cal reactions [^0]. Processes such as A + 5 ^ 0, where diffusing chemicals react irreversibly 
are believed to have many applications in physical, chemical, and biological systems. Par- 
ticular attention has been paid to cases where a reaction front is formed between regions 
dominated by A or i? particles. Such a situation can arise in the case where the two species 
are initially entirely segregated |3HlO| , or alternatively, and more simply, in a steady state sit- 
uation, where A and B particles are injected at equal rates at opposite boundaries |l9|JTT|-[l3 



In this letter we study the latter model, in the case of equal diffusion constants D for the 
two species. The simplest description of these systems is provided by the inhomogeneous 
mean-field rate equations for the particle densities a and b, where it is assumed that the 
reaction rate R = \ab: 



da 

at 



= DV^a - Xab (1) 



II = DV% - Xab. (2) 

These equations predict a reaction front width which scales as Wmf ~ {\J/ D"^)^^^^, where 
J are the (equal) imposed currents of A and B particles at the boundaries. However, it 
is well known that below a critical spatial dimension dc = 2 P JT^ , p!3| microscopic density 
fluctuations become relevant, and as a result the mean field approach breaks down. For d = 1 



Cornell and Droz |12| have suggested that the fluctuations modify the scaling of the width 



to w ~ {J / D) ^1"^. Numerical simulations ||TO|JT^ have broadly confirmed these conclusions. 



However, there has been recent controversy in the time dependent version of the model, with 
initially separated reactants, over the existence of multiscaling in the spatial moments of 
the one dimensional reaction front @,[10|- More recent simulations by Cornell [jlO[ have also 
indicated that the one dimensional profile is accurately described by a Gaussian. Hitherto 
this result has not been understood. 

Theoretical approaches to understanding the crucial role played by fluctuations have 
centered on mappings of the microscopic dynamics of the react ion- diffusion system onto 
a quantum field theory P JT3| , p!^ . This has allowed the effects of fluctuations to be sys- 



tematically included by summing sets of Feynman diagrams. Renormalization group (RG) 
techniques have then been employed to form a perturbation expansion in e — 2 — d. These 
calculations have confirmed the modified scaling in o? = 1, as well as pointing to the existence 
of power law tails, both in the densities and in the reaction front: 

R = AD{J/D)^\x\^^-^''>/^^ exp [-B{J/D)^x\^^-''>/^) 

+CD^J-'\x\-'+^' + ... (3) 

Here A, B, C are universal dimension dependent constants. One consequence of the power 
law tails is that sufficiently high order spatial moments of an evolving time dependent 
reaction front (with initially segregated reactants) should exhibit multiscaling. These re- 
sults were derived on the basis of a finite reaction rate, which under the RG was found 
to flow to a universal 0(e) fixed point. However, previous simulations of one dimensional 
reaction-diffusion systems have employed an infinite reaction rate, which enforces complete 
segregation between the two species. Furthermore only single occupancy of a given lattice 
site has previously been permitted. In this letter, we relax these restrictions by simulating 
a system with both multiple site occupancy and an adjustable, finite, reaction rate A. This 
model is closer to that used in the analytic RG calculations. 

Our model consists of a one dimensional lattice with L sites, on which particles of types 
A and B are located. In addition to this the model features reservoirs containing either A or 
B particles. The total number of particles of each species was set equal: Na — Nb — N/2. 
Three distinct processes take place: 

(1) A and B particles located on the lattice hop to neighboring sites in each direction 
with a hopping rate h, which we set equal to 1 (corresponding to = 1 in the continuum 
theory) . 

(2) Each A particle can react with each B particle on the same lattice site, with a reaction 
rate A. After each reaction both particles are removed from the lattice and placed in their 

respective reservoirs. 

(3) Each A {B) particle in the reservoir is inserted onto the leftmost (rightmost) site in 

3 



the lattice with an insertion rate i. Clearly J = Nres-h where Nres is the number of particles 
in the reservoir. The purpose of these reservoirs is to break up correlations between particle 
annihilation inside the reaction front and particle reinsertion at the boundaries - the larger 
the reservoirs the smaller the correlations. The same effect can be achieved by increasing 
the system size, but this is computationally far less efficient, as in that case particles have 
to hop large distances before they can reach the reaction zone. 

We carry out the simulations with rare-event dynamics (RED). In this approach no fixed 
time increment is present. First, in a specific configuration, a list is made of all the distinct 
events that might change the state of the system: 4L — 2 events for A and B particle 
hops, L events for recombination of a pair, and 2 events for insertion of a particle of either 
type; altogether 5L events. For each event ej, a rate rj is calculated. Each step in the 
RED simulation now consists of incrementing the time scale with At = 1/ J2j{'^j), and then 
allowing selection (and execution) of an event. The probability that event Cj is chosen is 
equal to Pi^n/ Y.j{rj)- 

For an efficient implementation, a binary tree of events is constructed, where each branch 
contains one event and has a weight equal to the rate of that event. The weight of a parent 
node is equal to the sum of the weights of its children. As the root node contains the sum 
of all rates, the time increment At is easily obtained. For the selection of a particular event 
Cj with rate r^, we start in the root node, descend to one of its children with a probability 
proportional to its weight, and iterate. The selected event is then executed and the tree is 
updated. 

The initial configuration for each simulation consisted of linear density profiles for the A 
and B particles, which decreased from the left and the right hand edges to the system center. 
In all our simulations, we chose L such that no A particle ever penetrated the B-rich region to 
within 10 sites of the lattice boundary, and vice versa. Correlation and thermalization times 
varied with J, N , A and L, and the tails of the reaction front required more thermalization 
than the middle section. The necessary thermalization time never exceeded 10^ events. To 
be safe, we thermalized our system in all runs over 10^ events. 
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We consider first the regime < 1, where the mean field reaction front 

width Wmf ~ {\J / D"^)^^/^ is much larger than the predicted fiuctuation modified width 
~ {J / D)^^/"^ . In this case we expect that the behavior of the system should be close to 
mean field. Solving equations (|l]) and for a and 6, we find that the mean field reaction 
front Rmf = ^ab has the form Rmf = J{XJ/ D'^y/^S{[\J/ V^^^'^x), where asymptotically 
(for {XJ/D^y/^\x\ > 1) we have 

S ~ i[\J/D'Y/'\x\)'Uxp (-^([AJ/DT/'|x|)i) (4) 

Hence the mean-field solution predicts that measured data for the reaction front should 
collapse if R/[J{\J/DY/^] is plotted function of {XJ/D'^Y'^x, as shown in figure 2. In 
this case the number of particles N (which varied from between 1300 to 12000) was tuned 
to obtain the desired J. 

The collapsed data is in good agreement with the mean field prediction, although there 
is a slight tendency for our simulation data to lie to the right of Rmf, for the largest values 
of {XJ/D'^y/^x. In this region, where the number of minority particles is small, we expect 
that noise from the reaction front will again become important, leading to a widening of the 
profile. Note that these simulation results were found not to depend on the inclusion of reser- 
voirs in our model, implying that the existence of correlations between particle annihilation 
and reinjection was unimportant in this case. 

In the limit XJ-^/'^D-^/'^ > 1 the mean field solution predicts that the reaction front 
will become increasingly narrow. However, the simulations disagree with this assertion - the 
reaction front keeps a finite width even if A is made very large. Our analysis, described below, 
distinguishes two components of this width: one is intrinsic, and the other is caused by the 
ability of the center of the front to wander. The intrinsic width is calculable using the RG 
approach already outlined, whereas the front wandering can be understood by considering 
the fiuctuations in the field ip = a — b, whose zero may be taken as defining the center of 
the front. Including the effects of reaction front noise (which is relevant in one dimension), 
the field theory leads to the following equation for ip: 
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^ = DV^ + r^. (5) 

Here 77 is the reaction front noise, satisfying (77) = and 

(r7(x, t)r/(x', t')) = 25{t - t')S{x - x')R, (6) 

where the reaction rate at the wandering front, with width Wg, has the form R = 
{J/wg)S{x/wg). It is important to reahze that, while ip is on average equal to (a) — (b), 
its fluctuations are not the same as those in the density difference. This arises from the 
non-trivial commutation properties of the operators within the field theory. More details 



on this point (within the context of an A + A ^ reaction) can be found in This 
fact accounts for the non-conservative nature of the noise in equation (^. We may now 
decompose ip into its mean field part together with higher order Fourier harmonics: 

^ = -{J/D)X + Y^Xnit) COS (^^!1±^^ (7) 

for — (L/2) < X < {L/2). These corrections are the most general possible which both 
couple to the noise (i.e. the harmonics have a non-zero amplitude at a; = 0), and which 
are appropriate for the non-conservative nature of the noise. Furthermore the densities on 
the boundaries are kept constant by these additional terms. We can now insert the above 
expression into the noisy diffusion equation and Fourier expand the reaction front noise 
(which is concentrated near x = 0). In the large time limit we find 



, , 2V2J /•* , 







L2 



{f - 1) 



dt', (8) 



where now 



m) = 0, mr,{t')) = 5{t-t'). (9) 
Clearly ip{x = 0,t) = I]Xn(^) is a Gaussian random variable with 

(V^(0,t)^) - (^(0,t))^ = E(Xn(t)x™(t)) (10) 

n,m 

Inf-l (11) 



in the large time limit, with c a constant, and where the upper limit is provided by the finite 
width Wg of the wandering reaction front. Assuming that the fluctuations of ip are small in 
comparison with the system size, then the gradient oi ip aX x = Q remains approximately 
equal to —{J/D). Hence to leading order we expect the position of the zero of the ip fleld 

to be a Gaussian random variable with width Wg, given by the recursive relation: 

^ 1 

\n{cL/wg) 



Wg 



(12) 



From our simulation data (in the limit AJ^^^^-D^^^^ ^ 1) we have plotted Rwg/J as a 
function of x/wg, where in ( |12|) we used c = 0.5 (see flgure 3). The collapsed data is well 
described by a normalized Gaussian, with width 1, in good agreement with our theory. This 
indicates that the higher order non-Gaussian corrections to the distribution of the zero of 
the ip fleld are indeed small. Note that the logarithmic factor in (|12]) is essential for a good 
flt to the data. 

We can clearly see from flgure 3 that the wandering Gaussian dominates over the intrinsic 
proflle to form the overwhelming component of the front. Notice also that we flnd a basic 
(J/D) -1/2 scaling, in agreement with earlier predictions [^,|T2|,|TE[ . Previously, however, only 



the intrinsic part of the front was being analyzed, whereas we have been studying the 
wandering piece. The scaling agreement is simply a consequence of dimensional analysis 
- any quantity with the dimensions of length, which is independent of A, must scale as 
(J/D)^^^^. We may also generalize our calculation to the time dependent case, where a 
reaction front is formed quasistatically between initially entirely segregated reactants. In 
this situation we flnd Wg ~ t"{lnty^'^, where a = 1/4. The presence of the logarithm may 



explain the slow convergence found in measurements of the exponent a [|T0[ . 

If we wish to study the intrinsic component of the front in the non-mean fleld limit 
{\J-^/^D-^l^ > 1), we must now flnd a way of suppressing the dominance of the wandering 
Gaussian part. One way in which this can be achieved is to measure the reaction rate Rr 
as a function of \x — the distance between successive reaction events at Xp and x. This 
enables the intrinsic proflle to be studied, as the front center has little time to move on such 



short time scales. These reaction events are effectively uncorrelated, so that the relative 
reaction rate Rr is given by 

Rr{x) = J R{x)P{xp)5{x — Xp — x)dxdxp, (13) 

where P{xp)dxp = R{xp)dxp/ J R{xp)dxp is the probability that the previous reaction oc- 
curred between Xp and Xp + dxp. In our simulations of Rr we kept the insertion rate i small, 
and hence many particles were present in the reservoirs. We were therefore able to effectively 
simulate a much bigger system, with a large number of particles. This ensured that correla- 
tions between annihilation and reinsertion, which would otherwise have modified the intrinsic 
profile, were kept to a minimum. Our simulation data is shown in figure 4, which shows a 
convincing data collapse of Rr{J / D)^^^"^ / J plotted against \x — Xp\{J / Dy/'^ . As we are now 
studying the intrinsic profile, this result finally confirms the scaling predictions of P, p!2| , p^ . 
For large values of \x — Xp\{J / Dy/'^ we also find a tail which is consistent with the RG im- 
proved tree level prediction log(i?) ~ —B{J / DY^'^\x\ (implying log(i?r) ~ —B{J/DY^'^\x\). 
Our data provides no clear indication of the power law tails predicted in |jIB[. However, 



such tails would be very hard to see in measurements of the relative reaction rate i?,., as the 
power law exponent would be large (7 — 2e + 0(e^)). Simulations have also been performed 
in the mean field regime. In this case Rr/[J{\J/D'^y/^] was found to collapse when plotted 
against small values of {\J / D'^Y^'^\x — Xp\. However, for larger values, the collapse no longer 
worked well, probably due to the increased importance of fluctuations in the asymptotic 
regime. 

In summary, our ability to adjust the reaction rate has enabled us to find two regimes for 
the A + B ^ % front in one dimension. For XJ^^I'^D^^I'^ <^ 1 mean field predictions work 
well, whereas for XJ^^^'^D^^^'^ ^ 1 the front is dominated by a Gaussian profile, a result 
of fluctuation induced wandering. Our theoretical prediction for this shape agrees well with 
simulations. Finally, we have succeeded in studying the intrinsic profile, where the shape 
and scaling properties match previous RG calculations. 
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FIGURES 

FIG. 1. In our model, A and B particles are located on a one-dimensional lattice. Three 

processes occur: particles can hop one lattice site to the left or right with rate h = 1; each pair of 
A and B particles at the same site can react with a rate A, after which they are moved to their 
respective reservoirs; and A {B) particles enter the lattice on the left (right) sides with an insertion 
rate i. 

FIG. 2. Collapsed data in the regime \J~^I'^D~^I'^ <C 1. SoHd line: mean field prediction. 
Squares: simulation results for runs over 10^ events and i = 1000, with D = 1, X = 0.001, 0.01, 0.1, 
and J = 0.1, 0.2, 0.5, 1.0. 

FIG. 3. Collapsed data in the regime XJ~^/^D~^/^ » 1. Solid line: normalized Gaussian; Data 
points: simulation results over 10^ events with D = 1, X = 1000, i = 1000, J = 0.1,0.2,0.5,1.0, 
and = 100 (0), 1000 (+). 

FIG. 4. Collapsed data for the relative reaction rate in the regime AJ~^/^D~^/^ ^ 1. Simula- 
tions were for 10^° events, with D = 1, N = 200, 1100, 10100. 
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